/*
 * gradientSaturation.cpp
 *
 *  Created on: 14/05/2009,2012
 *      Author: rogerio
 */

#include "EBFV1_hyperbolic.h"

namespace PRS{

	double EBFV1_hyperbolic::calculateSaturationGradient(TMesh* pTMesh_data){
		double start = MPI_Wtime();
		int dim = pTMesh_data->getMeshDim();
		int ndom = pTMesh_data->getNumDomains();
		int i, j, k, nrows, ncols, dom;
		for (dom=0; dom<ndom; dom++){
			// reset gradient data grad_Sw = 0 for domain 'dom'
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			for (j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(j,0);
				for(i=0;i<dim;i++){
					vdata->grad_Sw[i] = .0;
				}
			}

			// calculate gradient: all edges
			pTMesh_data->edge_list[dom].getSize(nrows,ncols);
			for (i=0; i<nrows; i++){
				EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);
				VertexData* vertex_0 = edata->vertex_0;
				VertexData* vertex_1 = edata->vertex_1;
				double val = 0.5*(vertex_0->Sw + vertex_1->Sw);
//				fid << vertex_0->ID << "\t"<< vertex_1->ID << "\t"<< vertex_0->Sw << "\t"<< vertex_1->Sw << "\t"<< edata->Cij[0] << "\t"<< edata->Cij[1] << "\n";
//				fid << vertex_0->grad_Sw[0] << "\t"<< vertex_0->grad_Sw[1] << "\t" << vertex_1->grad_Sw[0] << "\t"<< vertex_1->grad_Sw[1] << "\t";
				for (k=0; k<dim; k++){
					vertex_0->grad_Sw[k] += val*edata->Cij[k];
					vertex_1->grad_Sw[k] += -val*edata->Cij[k];
				}
//				fid << vertex_0->grad_Sw[0] << "\t"<< vertex_0->grad_Sw[1] << "\t" << vertex_1->grad_Sw[0] << "\t"<< vertex_1->grad_Sw[1] << "\n";
			}
//			fid << "\ncalculate gradient: only boundary elements (2-D: edges or 3-D: triangles\n\n";

			// calculate gradient: only boundary elements (2-D: edges or 3-D: triangles
			pTMesh_data->bdry_elem_list[dom].getSize(nrows,ncols);
			for (i=0; i<nrows; i++){
				BdryElementData* edata = pTMesh_data->bdry_elem_list[dom].getValue(i,0);
				VertexData* vertex_0 = edata->vertex_0;
				VertexData* vertex_1 = edata->vertex_1;

//				fid << vertex_0->ID << "\t"<< vertex_1->ID << "\t"<< vertex_0->Sw << "\t"<< vertex_1->Sw << "\t"<< edata->Dij[0] << "\t"<< edata->Dij[1] << "\n";
//				fid << vertex_0->grad_Sw[0] << "\t"<< vertex_0->grad_Sw[1] << "\t" << vertex_1->grad_Sw[0] << "\t"<< vertex_1->grad_Sw[1] << "\t";
				for (k=0; k<dim; k++){
					vertex_0->grad_Sw[k] += ((5.*vertex_0->Sw + vertex_1->Sw)/6.0)*edata->Dij[k];
					vertex_1->grad_Sw[k] += ((vertex_0->Sw + 5.*vertex_1->Sw)/6.0)*edata->Dij[k];
				}
//				fid << vertex_0->grad_Sw[0] << "\t"<< vertex_0->grad_Sw[1] << "\t" << vertex_1->grad_Sw[0] << "\t"<< vertex_1->grad_Sw[1] << "\n";
			}

//			fid << "\ncalculate gradient: weighting by node volume\n\n";

			// calculate gradient: weighting by node volume
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			for (i=0; i<nrows; i++){
				VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(i,0);
				for (k=0; k<dim; k++) {
					vdata->grad_Sw[k] /= vdata->volume;
				}
//				fid << vdata->grad_Sw[0] << "\t"<< vdata->grad_Sw[1] << "\n";
			}
//			fid.close();
//			STOP();

			// compute saturation gradient on partition boundary nodes
			// pMData->unifyVectorsOnMeshNodes(pStruct->SwPPData->get_Sw_Grad,pStruct->SwPPData->set_Sw_Grad,pGCData,3300);

			// project gradients:
			ofstream fid;
			//fid.open("grad.txt");
			pTMesh_data->bdry_elem_list[dom].getSize(nrows,ncols);
			double innerp;
			for (j=0; j<nrows; j++){
				BdryElementData* edata = pTMesh_data->bdry_elem_list[dom].getValue(j,0);
				VertexData* vertex_0 = edata->vertex_0;
				if ( !pSimPar->isInjectionWell( vertex_0->flag ) ){						// Excluding injection wells
					innerp = dot(vertex_0->grad_Sw,edata->versor,dim);			// Parallel Component of The Saturation Gradient
					for (i=0; i<2; i++){
						vertex_0->grad_Sw[i] = innerp*edata->versor[i];
					}
				}
				VertexData* vertex_1 = edata->vertex_1;
				if ( !pSimPar->isInjectionWell( vertex_1->flag ) ){						// Excluding injection wells
					innerp = dot(vertex_1->grad_Sw,edata->versor,dim);			// Parallel Component of The Saturation Gradient
					for (int i=0; i<2; i++){
						vertex_1->grad_Sw[i] = innerp*edata->versor[i];
					}
				}
				fid << setprecision(8) << fixed;
				fid << vertex_0->flag << "\t"<< vertex_1->flag << "\t" << vertex_0->ID << "\t" << vertex_1->ID << "\t" << vertex_0->grad_Sw[0] << "\t" << vertex_0->grad_Sw[1] << "\t" << vertex_1->grad_Sw[0] << "\t" << vertex_1->grad_Sw[1] << endl;
			}
			fid.close();
			//STOP();
		}
		double end = MPI_Wtime();
		return end-start;
	}
}
